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Abstract 

A general approach to a solution of few- and many-body scattering problems based on a 
continuum-discretization procedure is described in detail. The complete discretization of 
continuous spectrum is realized using stationary wave packets which are the normalized 
states constructed from exact non-normalized continuum states. Projecting the wave 
functions and all scattering operators like t-matrix, resolvent, etc. on such a wave-packet 
basis results in a formulation of quantum scattering problem entirely in terms of dis¬ 
crete elements and linear equations with regular matrices. It is demonstrated that there 
is a close relation between the above stationary wave packets and pseudostates which 
are employed often to approximate the scattering states with a finite L 2 basis. Such a 
fully discrete treatment of complicated few- and many-body scattering problems leads 
to significant simplification of their practical solution. Also we get finite-dimensional 
approximations for complicated operators like effective interactions between composite 
particles constructed via the Feshbach-type projection formalism. As illustrations to this 
general approach we consider several important particular problems including multichan¬ 
nel scattering and scattering in the three-nucleon system within the Faddeev framework. 

Keywords: quantum scattering theory, discretization of the continuum, Faddeev 
equations, multichannel scattering 
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1. Introduction 

The formulation and practical solution for scattering problems in few- and many-body 
systems entails on many difficulties caused by several reasons: (i) complicated boundary 
conditions, especially above the breakup thresholds; (ii) separation of the proper physical 
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solution from unphysical (e.g. forbidden by symmetry requirements etc.) ones; (iii) high 
dimension of resulting integral or integro-differential scattering equations; (iv) presence 
of complicated moving singularities in the kernels of few-body scattering equations. 

All these difficulties have resulted in the current situation in the field when we are able 
to treat rather accurately many-body bound states in very complicated systems consisting 
of tens or hundred particles while we are still unable to treat accurately even three- 
or four-body systems with Coulomb interactions, the scattering of particles with inner 
degrees of freedom, e.g. the collision of a few molecules in quantum theory of chemical 
reactions, or a triple collisions of three- or four nuclei etc. Moreover, even a solution 
of conventional three-nucleon scattering problem with realistic interactions requires an 
extensive usage of supercomputers [i[. It hampers strongly the further progress in this 
important area. Thus, one can formulate an important goal: not only to simplify greatly 
the solution of some few- and many-body scattering problems, but also to develop a new 
method in scattering theory which would be quite in a spirit of the bound-state problem 
solving. So that the whole treatment of few- and many-body scattering would be close 
to the routine bound-state calculations in quantum chemistry or in nuclear physics. 

On this way, many L 2 -type methods have been suggested and widely developed in the 
past and even now they are actively progressing. Among these the following methods can 
be mentioned: the R-matrix type approaches 0, the equivalent quadrature technique 
d, Q and the Stieltjes-Tchebyshev moment-theory approach [1,0, the Lorentz integral 
transform [ 1 ], the harmonic oscillator representation [ 1 ] and the J-matrix approach 0 , 
the continuum-discretized coupled-channel method (CDCC) Einiiii, the convergent 
close coupling approach [ll| , different realizations of the complex scaling method and 
others (see also the recent review ll]). All these approaches are utilized rather actively 
in atomic, chemical and nuclear physics. 

In this paper we suggest another general approach for solving a few-body scattering 
problem based on the usage of special normalized states, the so-called stationary wave 
packets (WPs) or eigendifferentials. The concept of eigendifferentials was introduced 
long ago by pioneers of quantum physics: H. Weyl [H |. E. Wigner 0 , H. Bethe [10 
and others to treat non-normalizable continuum states (which do not belong to a Hilbert 
space) in a framework of the standard theory of Hermitian operators in a Hilbert space. 
The idea of stationary wave packets turned out to be very fruitful at early stages of the 
quantum theory development (see e.g. the detailed discussion in the classical textbooks 

lua). 


In our previous papers [U, [^, [H, 24, 25, 26, 27, 28, [10 we proposed to use such 
normalized WP states as a very convenient discrete basis, onto which one can project 
out all the scattering operators and wave functions. Because WPs are related explicitly 
to exact scattering wave functions, many complicated operators such as free or channel 
resolvents have analytical finite-dimensional approximations in the WP representation. 
This leads immediately to a replacement of the scattering operators by finite matrices, 
so that a solution of a few-body scattering problem is reduced to solving simple matrix 
equations with regular matrix elements instead of multi-dimensional integral equations 
with singular kernels in a conventional formulation. In addition, our approach results in 
a very effective numerical scheme for a practical solution. Thus, in the present paper 
we report the main results of the general WP approach and describe how to reformulate 
complicated singular few-body scattering equations into a form which allows to solve 
them straightforwardly on a serial personal computer in a rather short time. 
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The work has the following structure. Section 2 is dedicated to the description of 
the eigendifferentials (stationary wave packets) and their interrelation to conventional 
non-normalized scattering wave functions, on the one hand, and discrete pseudostates 
resulting from a standard diagonalization procedure for the Hamiltonian matrix, on the 
other. Section 3 includes a detailed description of different properties and features of 
the stationary wave packets. In Section 4 we present a discrete version of the scattering 
theory based on the WP representation for the scattering operators and wave functions. 
Section 5 is devoted to the multichannel scattering while the next Section 6 contains a 
general description of the three-body wave-packet formalism. In Section 7 a solution of 
the composite particle scattering off a target-nucleus is discussed. A formulation of the 
Faddeev equations in a discrete form and their practical solution are given in Section 8. 
We summarize our main results in the last Section 9. For the reader’s convenience we 
added the Appendix which contains explicit formulas for the channel resolvent eigenvalues 
in the WP representation. 

2. Discrete representation for the scattering states 

Let us introduce some notations. We start here with two-body scattering problem in 
a stationary formulation. Assume that the potential v is local and spherically symmet¬ 
rical, so that the angular momentum is conserved and the total Schroedinger equation 
is reduced to radial equations for a fixed value of the angular momentum 1. Below we 
will deal mainly with the radial parts of the wave functions which justify the partial 
Schroedinger equation; 

+ v{r)^'tp'-{r,E) = E'ip'^{r,E), (1) 

where we use units with h = 1 and ^ is the reduced mass. We will omit the index I 
wherever possible until Section 6 and also other indices related to spin, isospin etc. until 
the cases they are to be detailed. 

The operator in the left-hand side of eq. o is the total Hamiltonian h, which can be 
written in the form: 

h = ho +v, (2) 

where ho is a free motion Hamiltonian. Further we will write all the relationships in 
an operator form using the Dirac notations, without specifying the representation. It is 
essential that Hamiltonians ho and h are Hermitian operators (in an appropriate Hilbert 
space) and have a simple continuous spectrum [0,oo). 

The total Hamiltonian h may also have a discrete spectrum describing bound states of 
the system. The corresponding bound-state wave functions of h are denoted as {|'0n)}^=i 
while the continuum functions are denoted as |'i/'(7f)). They satisfy the Schroedinger 
equation (an operator form of the eq. ©) 

h\i,{E)) = E\i,{E)) (3) 

and the usual orthogonality condition 

{i,[EmE'))=5[E-E'). 
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<p + 


dr^ ' J.2 


( 4 ) 





In our study, we also need the scattering wave functions j'i/'g) = 
to the delta-function of momentum q = 

(V'-jIV'?') = S{q-q'). 


normalized 


( 5 ) 


2.1. Description of a continuous spectrum in terms of eigendifferentials 

A description of continuum meets some difficulty in a stationary formulation of a 
quantum mechanical problem. As is well known, continuum wave functions for the 
total and free Hamiltonians, h and /iq, do not have a finite normalization unlike to 
bound-state wave functions and thus they do not belong to a Hilbert space. Strictly 
speaking, such non-normalizable states are not eigenfunctions of h and ho and some 
basic properties of Hamiltonians (the Hermiticity property) or even an orthogonality 
of continuum wave functions could not be proven in a normal manner with such non- 
normalizable set 17, 1^, 2^. 


Nowadays this formal problem is avoided by using the Dirac delta-function (|1|) and 
by generalization of a Hilbert space to a rigged Hilbert space. However, in the past, 
physicists often employed another formal trick (see e.g. [171, [l^ [lOl l20|) based on the 
Weyl’s eigendifferential concept developed when studying the spectral theory for singular 
differential operators 


161 . 


The idea is to introduce an interval with a small width AE for each value E from 
the continuum and then to construct the so-called eigendifferential, i.e. integral of the 
continuum wave function |^(£')) over the interval: 


mE,AE)) = 


f'E-\-AE 


dE\ij{E)). 


( 6 ) 


It is easy to see that such state is normalizable (due to an integration) and belongs to 
a Hilbert space. In this way, one can generalize the conventional definition of a state 
normalization: the state \ip{E)) is treated as normalizable if its eigendifferential has a 
finite norm 


181. 


To treat a whole continuum, it should be divided into non-overlapping intervals (i.e . 
it should be discretized). The system of eigendifferentials forms the orthonormal set [20| : 


{^{E',AE')\i:{E,AE)) = 


AE, for the same intervals, 

0, for different intervals 


( 7 ) 


It is important to stress that the set of such eigendifferentials give a diagonal represen¬ 
tation for the Hamiltonian: 


{i;{E\AE')\h\i;{E,AE)) 


E + \AE, for the same intervals, , , 

Q frM' rliflPo-rcinf i n^-d-rTT-o 1 o ' ' 


AE 1 0, for different intervals 

Besides this, they have a finite overlap with initial (non-normalized) states 

I, E' €{E,E + AE), 


{i;{E'mE,AE)) = 


0, E' i{E,E + AE) 


( 9 ) 


A complete system of eigenfunctions for the Hermitian Hamiltonian h in the eigendif¬ 
ferential formalism is built from bound-state wave functions and eigendifferentials 
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Elis- In other words, for an arbitrary wave function |$) one can write down the 
following expansion; 


Nb 

|$) = ^C„|V'^) + ^C(i?)|V^(i?,Ai?)), (10) 

n—1 

where the sum over an infinite but accountable set of eigendifferentials has rather a 
symbolic meaning. However further one can pass to a limit AE —> 0 which leads imme¬ 
diately to a conventional expansion of the arbitrary function |$) over the bound-states 
and functions of the continuum: 


/.OO 

1$) = ^ + / dEC{E)\i;{E)) (11) 

n=l “'O 

The discrete sum in eq. ra goes to the integral expansion in eq. (Hill over continuum 
wave functions \’tp{E)) in the sense of a Riemann-Stieltjes integral [20[. Also it can be 
shown pojl that the finite normalization condition for eigendifferentials o leads just to 
a delta-function normalization of continuum states ®. 

Today such a way of definition of the Hamiltonian eigenfunctions for a continuous 
spectrum practically is not used, as was noted above. But in our approach we use such 
eigendifferentials as basis states to solve a scattering problem in a discrete representation. 
Below we demonstrate that this “classical” way is turned out to be extremely fruitful for 
practical calculations in few-body scattering problems. 

There is also the another very popular concept for a treatment of a discretized con¬ 
tinuum: the pseudostates. It will be demonstrated below that these two ways of dis¬ 
cretization: via eigendifferentials and via pseudostates are closely interrelated. 

2.2. Pseudostates and L 2 discretization of a continuum 

As is known, pseudostates arise in an approximate treatment of a continuous spectrum 
of a quantum system based on a projection of scattering wave functions into an orthogonal 
L 2 basis set {\4>n)}n=i of ^ finite dimension N. 

If one substitutes the expansion of a Hamiltonian eigenfunction in the basis functions 
1^) = into the Schroedinger equation ([3|) one gets the following system of 

algebraic linear equations for the expansion coefficients {Cn}n=i 

N 

Y^[hrbk-E6nk]Ck=0, n = l,...,N, (12) 

fc=i 

which corresponds to the eigenstate problem for the Hamiltonian matrix with the ele¬ 
ments hnk = {4’n\h\4>k)■ After a diagonalization of the Hamiltonian matrix one gets a set 
of eigenvaluelil and eigenfunctions \ifn) = Y^k=i^k\^k) {n = 1,...,N). This 

discrete set can be divided on two groups: the eigenfunctions with lower eigenvalues 
describe the system bound states if they exist (in the case of a short-range potential 


^ Below we denote eigenvalues by an upper asterisk index to distinguish them from discretization 
interval end-points. 
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their number is finite) while the remained eigenstates correspond to the continuum wave 
functions at discrete energies e*. Unlike to the exact continuum functions these discrete 
eigenstates have a finite norm and their eigenvalues define the discretized continuum. 
That is why they are usually referred to as pseudostates of the Hamiltonian h. 

For the above pseudostates, the conditions similar to those given in eqs. (HHHl) for 
eigendifferentials are valid: 

(V’nIV'fe) “ Snkt 

{i>n\h\llJk) = in^nk ^ 

An employment of pseudostates instead of exact scattering functions in practical scatter¬ 
ing calculations meets the well known difficulty related to the difference in their normal¬ 
ization. The functions of pseudostates behave in some inner region like the scattering 
wave functions but, in sharp contrast with the latter, they vanish in the asymptotic 
region. Nevertheless, the functions of pseudostates together with bound-state wave func¬ 
tions form a complete set of basis functions, and this set can be used to approximate the 
spectral expansions of the scattering operators or scattering wave functions. 

2.3. An equivalent quadrature concept 

Very often it is necessary to find not operator itself, but only its matrix elements 
with normalized functions |$). The example for such an employment is the calculation 
of transition probabilities between the bound states of subsystem via the intermediate 
continuum states, in particular, the calculation of response functions in the Lorentz 
integral transform (LIT) method Q or the calculation of f-matrix elements on the basis 
of spectral expansion of the total resolvent of the system 

Let us consider a matrix element of some operator B = f{h) which is a function of 
the Hamiltonian h using the spectral expansion of the operator H in a complete set of 
the Hamiltonian eigenfunctions: 


Nb 

($| H | cI >) = ^/( U „) K < 1>|^(()|2 + 

n—1 

As is evident, the spectral expansion da includes not wave functions themselves but 
their positive quadratic forms \{^\tfj{E))\‘^. So that, for such quantities there are no 
problems related to the normalization of scattering wave functions and they are finite 
functions of energy. Further, the integral over continuous spectrum can be approximated 
by a discrete sum over respective pseudostates: 


pOO 

/ dEf{E)mi;{E))\\ (14) 

^0 


Nt N 

(cI>|H|ci>)«^/(u:)|(ci>|^0)|2+ ^ /(F;:)|(cI>|V;„)p. (15) 

n—1 n—Nb-\-l 

Thus, such a “spectral expansion” over pseudostates can be treated as some effective 
quadrature for the integral over continuum in eq. (O- This approach has been sug¬ 
gested earlier for solution of different particular problems in atomic physics mm 
related to finding matrix elements of Green functions. E.g., a conception of an equivalent 
quadrature (EQ) was introduced within the framework of the moment-theory approach 
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for calculation of non-negative spectral densities ill- The mesh points of the EQ coin¬ 
cide with pseudostate energies and corresponding weights are determined formally from 
relations between discrete and continuous quantities: 

n = Nt + l,...,N (16) 

which are valid at discrete pseudostate energies e*. 

One important result of the equivalent quadrature technique is the fact the quadrature 
weights do not depend upon the function |<I)) and are determined actually by pseudostate 
wave functions and by the spectral density of the particular system. So that, the values a;„ 
are just the coefficients for transitions from normalizable states \tjjn) to non-normalizable 
functions |^(e)))). 

In general, finding the weights of the equivalent quadrature is highly non-trivial prob¬ 
lem because the exact scattering wave functions \tjj{E)) are still unknown. However, for 
many problems when one needs only complete spectral sum, knowledge of these weights 
is not required. For example, the pseudostate functions are used for inclusion of interme¬ 
diate continuum states for excitations of projectile or target in description of composite 
particle scattering in atomic [l^ and nuclear [l^ physics. In these cases, as a rule, the 
pseudostates are employed in expansion of a total wave function of the system to perform 
a coupled channel reduction of the initial many-body problem. 

However, it is still possible to find scattering observables just from individual pseu¬ 
dostate wave functions if the EQ is known. In fact, if to write down an expression for 
the elastic scattering amplitude: 

A{E) = (MEMm)), (17) 

it is easy to see that the function v\'tjjo{E)), where |'0o(7?)) is a free motion function, 
is square-integrable in case of a short-range potential v. So that, one can employ the 
properly normalized pseudostate instead of the exact scattering wave function: 

An « (V'0«)k|V'n). (18) 

To find direct relation between “continuous” A{E) and “discrete” amplitudes, one 
can use the EQ definition which results in: 

A{En) ^ (19) 

y/^n 

Actually, the problem of evaluation of is related to finding a spectral density in the 
continuous spectrum of h. One of the methods is to evaluate moments of the energy 
distribution using spectral expansion in pseudostates. Then, the equivalent quadrature 
weights can be found from a solution of complicated set of nonlinear algebraic equations 

However, to find the EQ, the eigendifferential formalism can be used. 


2.4- Pseudostates as approximations for eigendifferentials 
In our previous works 


22, 23 


the following way for finding the coefficients of 
transformation from normalizable functions to non-normalizable ones has been suggested. 
This method is based on the observation that pseudostate wave functions are very similar 
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to eigendifferentials. In particular, both types of functions behave in the inner region like 
an exact scattering wave function, whereas in the asymptotics they vanish differently in 
accordance with their own properties. Moreover, functions from both sets are normalized 
and the Hamiltonian matrices in both sets are diagonal. 

It means one can treat pseudostates as approximations just for eigendifferentials 
rather than for exact non-normalizable scattering wave functiont^. Of course, this ap¬ 
proximation is valid only when eigenenergies of both states m and (ESI) are coincided. 

In Fig [1] we compare pseudostates of the free Hamiltonian found in the Gaussian basis 
(see details in ref. [2g) with exact eigendifferentials (normalized to unity) constructed 
for the same distribution of eigenenergies e*. It is clear that both type functions are 
nearly coincide in some inner domain. 





Figure 1: Comparison of pseudostates (solid curves) of the free motion NN Hamiltonian obtained in a 
basis of Gaussians with exact eigendifferentials (dashed curves) constructed for the same distribution of 
eigenenergies of the free Hamiltonian discretized continuum. 


Further, it is easy to see that the coefficients of transformation to non-normalizable 
states (see the definition Q) are given simply by the widths AE in the discretization 
procedure for eigendifferentials: 


_ {ME)\v\'fp{E,AE) 


( 20 ) 


Using the direct correspondence between eigendifferentials (with the normalization fac- 


^It is worth noting that the interrelation of pseudostates and eigendifferentials has been also proved in 
refs. Em on the basis of the CDCC method. In this approach these two types of states are considered 
to treat a composite particle scattering off heavy target. However, direct relations between them have 
been established for the equidistant discretization distribution only. 


9 








tor) and pseudostates 


yj AEn 


(V'0(4)k|^n), 


and comparing the eqs. m and (uni, one gets an approximate relation: 


LOn ~ 


( 21 ) 


( 22 ) 


So that, we find an amazing result: the coefficients of transformation from pseudostates 
to exact continuum wave functions are determined, first of all, by the distribution of 
discrete pseudostates energies e* rather than by the form of basis functions |(()„) them¬ 
selves. Actually the form of basis functions determines implicitly the pseudostate energy 
distribution: it will be different for various basis sets. 

It would be very appropriate to give here some clear example. Let us consider a free 
particle motion in a spherical box with radius R. The solution of a radial Schroedinger 
equation for S-wave in this case is an infinite set of wave functions (normalized to unity) 

V' 0 n(?’) = Vising*r, = ^,n = 0 ,..., ( 23 ) 

instead of exact free motion states: 

i’Oqir) = \ -sinqr. ( 24 ) 

V TT 

For this case the weights of the EQ do not depend on the index n and can be found 
explicitly: 

OJn = ( 25 ) 

It is clear that this quantities coincide exactly to the discrete momentum differences 
dn = qn+i — <ln (momentum analog of the energy widths in ([51)). In other words, one can 
determine the weights of the EQ directly via the discretization width 

UJn = dn ( 26 ) 

rather that via the size of the box R (although they are interrelated). 

When we are dealing with an L 2 -discretization using an arbitrary basis, this can be 
treated as a solution of the scattering problem in a box with blurred boundaries but these 
boundaries are difficult to determine explicitly. However they do impact on a discrete en¬ 
ergy distribution which is used for finding approximate weights or widths of discretization 
intervals. So that, eventually in order to derive the values of scattering observables from 
the discrete energy distribution one should determine the spectral partition parameters 
using eigenvalues e* of pseudostates obtained from the Hamiltonian matrix diagonaliza- 
tion. In general case, the energy interval widths AE„ would be different for different 
n. 

Thus, our main idea here is that pseudostates are approximations for exact eigendif- 
ferentials. Surely, the above argumentation is not strict but plausible. It is also necessary 
to have similarity between pseudostates and eigendifferentials at least within the inter¬ 
action range. It means that the size of the pseudostate basis should be rather large and 
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the interval widths to be rather small. Unfortunately, it is difficult to give more strict 
conditions for such a replacement here. Some estimates can be found in ref. [^ . 

Below we formulate a fully discrete version of a stationary scattering theory using the 
eigendifferential (wave-packet) formalism, while in practical realizations we replace often 
exact eigendifferentials with pseudostates found from the Hamiltonian matrix diagonal- 
ization on an appropriate basis. Both these points make it possible to simplify greatly 
the practical solving of few-body scattering problems. 

3. Stationary wave packets and their properties 

In the approach which will be used everywhere in this work, the eigendifferentials 
will be employed as primary basis functions upon which all wave functions and opera¬ 
tors will be expanded. In other words, we will treat scattering problems in a discrete 
representation. 

We rename eigendifferentials as stationary wave packets because the name seems to 
us more physical and clear for understanding. The approach developed — we refer to it 
as the wave-packet continuum discretization method (WPCD) — takes the advantages 
of the above L 2 techniques and simultaneously gives an accurate solution of few-body 
equations. 

3.1. Basic definition of stationary wave packets 

It is important to stress that the stationary wave packets are constructed from exact 
wave functions of continuous spectrum as their eigendifferentials (normalized to unity). 
But WPs are considered here as ordinary functions and thus the details of a partition 
procedure for a continuous spectrum have to be defined. 

Let us restrict the continuous spectrum of the free Hamiltonian hp [0, 00 ) with some 
maximal value of energy Umax and then divide the whole interval [0, Umax] into finite 
number of non-overlapping intervals (where £q = 0 and £n = Umax) which 

we recall, according to ref. [l^, as discretization bins. We assume here that the value 
Umax is finite although sufficiently large in order to provide convergence to the exact 
solution of the problem to be treated. Every such an energy bin corresponds to the 
interval [g^-i, qi] on the momentum axis q where q = y/2pE and /i is the reduced mass. 
For the sake of convenience we will denote both the momentum and the energy bins by 
a symbol 'Si {i is an interval number). Also we denote widths of momentum and energy 
intervals as follows: 

di — qi — 1 , B)i = £i Si — i. 

Further we consider the complete set of continuous spectrum states {ipoq) for the free 
Hamiltonian ho (for each partial wave 1) normalized as in eq. (O. Now we define the free 
stationary wave packets as integrals of free motion wave functions over the discretization 
bins as in eq. (IH): 

= d,qf{q)\-ipoq), i = l,...,N, (28) 

V-Di Jx>i 

where f{q) is some weight function and Bi is a normalization coefficient related to each 
other by the formula: 



( 29 ) 




By modifying the weight function f{q) one can get different types of WPs. In particular, 
the case f{q) = \Jq/^ corresponds to the energy wave-packets which coincide with 
eigendifferentials (jH]) up to normalization factors. The case f{q) = 1 leads to momentum 
WPs etc. However it is important to note that when the partition used is sufficiently 
dense (i.e. for small bin widths) the particular form of the weight function does not 
matter. Moreover in practical calculations one can replace energy WPs with momentum 
ones (for the same partition) and errors related to this replacement are very small. 

According to the definition (E51) it is easy to show that the states |xi)A ^ form an 
orthonormalized set (^ : 


{xi\xj) = 6ij, = 


(30) 


Then it is straightforward to get a formula for the overlap integral between WP states 
and initial plane waves: 


/(?) 


9 e S, 


0, q> qN- 


(31) 


The free Hamiltonian matrix (and also the matrix of any operator depending upon 
or just commuting with it) is diagonal in such a basis while the respective eigenvalues 
depend on the weight function: 


{xi\ho\xj) = 4^ / |/(g)p|-dg. 




i J'Si 2/r 


In the case of energy WPs, the eigenvalues of ho are just the bin midpoints: 

{xi\ho\xi) = £* = +£i), 

and in the case of momentum WPs, the eigenvalues are: 

+ 9?)- 


(32) 


(33) 


(34) 


Such WP states will be used to carry out all our calculations. From practical point of 
view the wave-packet set can be used in any scattering calculations on the same footing 
as usual discrete L 2 bases like the harmonic oscillator basis or a basis of Gaussians. 
However, in contrast to those, such WP basis states have a direct relation ra to exact 
continuum wave functions. 


3.2. Behavior of 'wave-packet functions in configuration and momentum spaces 

The WP functions have an interesting behavior both in coordinate and momentum 
spaces. Consider for simplicity the case of f{q) = 1. It is easy to show 22, ^ that in 
coordinate space WP functions behave as follows (here we use the reduced radial part of 
a wave function for s-wave): 


Xi{r) 



sin(g* r) 
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Figure 2: Coordinate behavior of s-wave WP states Xi{r) for different values of bin widths in ratio to 
momentum value q\ di/q = 0.25 (a), di/q = 0.1 (b) and di/q = 0.05 (c). 


where q* = + qi) is a midpoint of the interval 2)^. So that, free WPs are nearly 

coincide with non-normalized continuum wave functions in the area r <C = ^. As is 
seen from eq. (123, the “packetting” procedure defined in eq. results in an additional 
decreasing factor which just provides the wave function with finite normalization. 

In Fig. [2] the coordinate behavior of three WP states (I35|) are compared for different 
ratios of bin width di and its midpoint value q*. It follows from the Figure that the radial 
behavior of the WPs resemble the characteristic picture of beats (i.e. dual-frequency os¬ 
cillations) and the smaller ratio di/q* corresponds the slower decay of the WP function. 
Therefore such a basis functions can be employed to approximate the non-normalized 
continuum wave functions on longer distances. This property is crucially important for 
few-body scattering calculations because the Faddeev components of the total wave func¬ 
tion have a specific long-range asymptotic behavior which must be reproduced correctly. 

In the momentum representation, the free WP functions are simple step-like functions: 
they do not vanish only for the momentum interval coinciding with their bin, i.e. only 
for the on-shell region. In this area they are fully determined by the weight function 

fid)-- 


i ^ f{q){S{q-q^-l)-^{q-q^)) 


(36) 


where the 6 is the Heaviside step-like functior^. When f{q) = 1 all the wave functions 
in the momentum WP representation take a histogram form. Being generalized onto a 


^Hear the ^-functions is defined as follows: 


(J > 0 
(J < 0 
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few-body case, the free few-body WP basis functions are built as a direct products of 
step-like functions along every independent momentum variable. In this way, the whole 
momentum space is replaced by a finite momentum lattice. In this sense we will refer to 
this basis as a lattice basis. 


3.3. Wave packets as a basis 

For further investigations, the important questions are: does the WP set form a basis 
and how one can prove a convergence of the results with increasing a basis dimension. It 
should be noted here that for the set of pseudostates the problem of completeness does 
not exist because the pseudostates jointly with the bound states form the same linear 
span as the initial L 2 basis set. So that, the set of pseudostates and bound states form 
a basis which becomes to be complete in the limit N ^ 00 . 

It is clear that, unlike to the eigendifferential case when AE —?► 0, the WP set does 
not form a complete basis in a full Hilbert space if bin widths are kept to be finite. 
But one can introduce the WP subspace Hp in which the orthonormal set {|a:i)}^i, of 
course, forms a basis. The projector onto the WP subspace is defined as usual: 

N 

P = XI 

i=l 

Thus below we will refer to the WP set as a WP basis keeping in mind it is the basis in 

'Up. 

One has to make some estimates concerning the convergence of the results obtained 
in "Hp to exact scattering-problem solutions with increasing the basis dimension. 

The exact unit operator /, which can be written in the following form: 


N 

I = X] 

i=l 





\'fp0q){'fp0q\dq, 


(38) 


is replaced in the WP approach with the wave-packet projection operator p. This re¬ 
placement implies two following approximations: 

(i) the infinite continuous spectrum is truncated with the maximal momentum value 
qn and the residual integral is neglected; 


(ii) exact partial spectral projectors li are replaced with the WP partial projectors 

\Xi){Xi\. 

Surely, these points are not valid in a full Hilbert space. But keeping in mind numerical 
applications of the method one can compare the mean values of operators I and p in 
some L 2 normalized state |$) having an effective range tq. To satisfy to the following 
conditions (corresponding to the above points (i) and (ii)) 

N N 

($|p|$)«($|X^*l<J>), (<i>|/.|<i>)«($lX^*l‘5’) (39) 

i=l i=l 

one has to have sufficiently small widths di and sufficiently high maximal momentum 
value qn (^ : 

di < — < giv, i = l,...,N. 
ro 
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(40) 


With these restrictions one can choose for the practical realization the proper condi¬ 
tions for the momentum bin partition Then one can check the convergence of 

the results with increasing a wave-packet basis dimension, when the bin widths become 
smaller and smaller, and the maximum value qn becomes higher and higher. This proce¬ 
dure simulates changes in eigenvalue {e^} distributions for pseudostates with increasing 
a basis dimension N. 

One can employ, for example, distributions (CTids) transforming a finite interval to 
whole numerical axis, e.g. the Tchebyshev grid j27j|l: 


= Qm tan 


2i- 1 


AN 


= 1 , 


,iV, 


(41) 


where is a scale parameter. 

Such grid points qi are convenient because one can study the convergence of the 
results with increasing the basis dimension as in the usual L 2 -basis case. One can expect 
that the result obtained in the WP representation approaches the exact solution with 
increasing the basis dimension N. 

3.4- Construction of WP states for the total Hamiltonian 

Similarly to the free Hamiltonian (ho) case, one can introduce a WP set for the 
total Hamiltonian h = v. We define a partition of the continuum for the total 

Hamiltonian into non-overlapping bins {Afc = with cq = 0 and cm = ^'max 

(or Afc = [qk-i^qk])- Let’s note that this partition might be different from the free WP 
partition {Si}. The scattering wave-packets are defined as integrals of exact scattering 
wave functions j'i/'g) over the intervals chosen: 

\zk) = [ dq w{q)\'iljq), (42) 

VL'fe JAk 

where Ck and 'w{q) are the normalization coefficients and weight functions respectively. 
They have the same properties as free WPs. However, the total Hamiltonian may have 
a discrete spectrum. In this case one has to add bound-state wave functions of h to the 
scattering WP set to obtain complete WP set corresponding to the total Hamiltonian h. 

In practical solution we have the basic problem how to construct these scattering 
WPs without solving the initial scattering problem. Here one can exploit an idea of 
similarity of eigendifferentials (WPs) and pseudostates. So that, one can use free WPs 
as a basis to construct scattering WPs for the total Hamiltonian h as its pseudostates. 
For this purpose we apply a diagonalization procedure for the total Hamiltonian matrix 
in a free WP basis. As a result we get a discrete sets of eigenvalues e*f. and respective 
eigenvectors \zk)- 

Since the free WP basis functions (in the momentum space) are step-like functions, 
the momentum dependence of all functions expressed via such a basis have a histogram¬ 
like form. An example of the momentum dependence for the bound state (deuteron) 
function in such step-like basis in comparison with the exact function for the Yamaguchi 
s-wave triplet NN potential is displayed in Fig. |31 

All the eigenstates corresponding to the area of continuous spectrum obtained in such 
a WP-representation can be treated as approximations for scattering wave packets \zk) 
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Figure 3: Comparison of the exact deuteron wave function in the momentum representation for the 
Yamaguchi potential (dashed curve) with its approximation in the lattice basis (solid line). 


defined above, provided the eigenvalues of Hamiltonian matrix coincide with eigenvalues 
of the exact wave-packets. This condition can be easily satisfied because one can construct 
a bin partition for exact scattering WPs “by hands”. 

As a result, we find very convenient discrete representation for the scattering WPs 
using the free wave packets: 


N 

\zk) ^ \zk) ='^Oki\x^). (43) 

i=l 

The validity of this approximation is illustrated by the same particular example of the 
Yamaguchi potential for which both exact solutions and also exact scattering WPs are 
constructed in an explicit form and thus one can compare them with approximations 
given by eq. (021) • In Fig- IHwe display the functions of two pseudostates in momentum 
representation (with indices k = 6 and 11) obtained in the free WP basis in comparison 
with the corresponding exact scattering wave packets. It is interesting that the exact 
scattering WPs (0^ (dashed curves in the Figure) are square-integrable functions, despite 
the fact that they have logarithmic singularities on the boundaries of the “on-shell” 
interval, i.e. the interval to which the WP energy belongs. It is clear from the 
above comparison that the pseudostates composed from step-like WPs reproduce quite 
reasonably the structure of exact scattering wave packets “on average”. 

Such a discrete approximation for the scattering WPs is extremely important for few- 
body scattering studies where one is able to build a few- and many-body WP basis not 
only for a free motion Hamiltonian but also for a few-body channel Hamiltonian. More¬ 
over, in this approach one can find a very convenient analytical form for the resolvent and 
off-shell t-matrix operators for a two-body subsystem by a straightforward Hamiltonian 
diagonalization procedure. 
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Figure 4: (Color online) The functions of pseudostates (fc = 6,11) obtained in the lattice basis (solid 
lines) in comparison with exact scattering wave-packets (dashed lines) for the NN spin-triplet Yamaguchi 
potential. 


Furthermore we will demonstrate henceforth how to generalize this convenient way 
to the general multichannel scattering. However, in general case of a multi-particle 
scattering one has to take into account a specific behavior of scattering wave functions 
in different asymptotic channels. In such cases we will formulate discrete versions for the 
Faddeev and many-body Lippmann-Schwinger integral equations, see Sections 6-8. 


3.5. Coulomb wave packets 

One of the big advantages of an employment of L 2 normalized WP states for a con¬ 
tinuum treatment is a possibility to consider long-range potentials like the Coulomb one 
quite similarly to short-range ones. It is because the essential singularities at small mo¬ 
mentum q peculiar to wave functions and transition operators for the Coulomb potential 
are averaged out and smoothed in the WP representation. 

As an illustration let us consider two-body scattering with repulsive Coulomb Hamil- 
toniaifl 

= + (44) 


where Zi and Z 2 are the particle charges, r is the distance between them. To treat the 
case of charged particles one can employ the Coulomb wave packets \xf) as the basis 
functions. Such Coulomb WPs may be built from the regular Coulomb wave functions 
F/(( 7 , r) by the integration over discretization bins 21| quite similarly to the general case: 


kf) = [ dqfiq)\Fi{q)). 

V JVi 


(45) 


^The Coulomb attraction can be also treated in the WP approach, however it needs in a separate 
study. 
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It can be shown straightforwardly that WP states (|45|) are normalized and thus one can 
construct them practically using pseudostates of the Hamiltonian (l4^ on some L 2 basis. 
Thus, in such a discrete representation one gets an interesting picture when one can 
expand the Coulomb WPs over finite set of free-motion WPs. 

This statement is illustrated in Fig.[5]where the exact Coulomb WPs for pp system are 
compared to respective Coulomb pseudostates found via diagonalization of the Coulomb 
Hamiltonian on the free WP basis, and free WPs themselves at the same energy. It is 
clearly seen from the Figure the Coulomb WP \xf) can be very accurately approximated 
by free WPs. 


x^ir) (fm-'“) 



Figure 5: The exact Coulomb WPs (dashed curves), the pseudostates found in the free WP basis (solid 
curves) and the free WPs at the same energy (dash-dotted curves) for pp system at three center of mass 
energies: Ec.m. = 0.03 MeV (a), Ec.m. = 0.133 MeV (b) and Ec.m. = 1.474 MeV (c). 


4. Discrete version of the scattering theory 


The explicit relation between continuum wave functions and their WPs makes it 
possible to develop a closed wave-packet formalism to treat scattering in a very convenient 
discrete representation. The complete formalism has been described in detail in our 
previous papers 26 | so that we present here for the reader convenience only the 


extract of the basic results. 


4-.1. The wave paeket space 

Let us call the linear shell spanned on basis of free WPs as an eigen wave- 

packet space of the Hamiltonian h^. The projection operator onto this space p is defined 
in eq. (1371) . 
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The property m allows one to find a finite-dimensional representation for any op¬ 
erator A = R{ho) which depends on Hamiltonian ho or just commute with it: 

N 

a = pHp = ^ (46) 

where corresponding eigenvalues Ai are given in an explicit form: 

E.g., discrete eigenvalues of the free Hamiltonian are given via eq. dsa). 

So that, using these properties one can build finite-dimensional analogs for differ¬ 
ent scattering operators in the free WP basis. In particular, the resolvent of the free 
Hamiltonian go{E) = [E + iO — ho]~^ is diagonal in the free WP representation: 

N 

qo{E) = pgo{E)p = ^ \xi)gi{E){xi\. (48) 

i=l 

Here the complex eigenvalues gi{E) are determined according to the general formula (H7)) . 

For the charged particle scattering, one may easily build the finite-dimensional ap¬ 
proximation for the Coulomb resolvent gc = [E + iO — hc]~^ in the Coulomb WP basis: 

N 

Sc{E)=J2\xf)g.{E){xf\ (49) 

i=l 

The Coulomb resolvent eigenvalues gi{E) are found from the same eq. (I47|) as for the 
free resolvent. The explicit formulas for these eigenvalues are given in the Appendix A. 

The eigenvalues gi{E) incorporate logarithmic singularities at the bin endpoints. To 
smooth such undesirable singularities and to convert our scheme to be completely dis¬ 
crete, one can employ an additional averaging procedure over the “on-shell” energy bin 
to represent the energy dependence of the resolvent. As a result, we get a purely finite¬ 
dimensional representation for the free resolvent which is free of any singularities at real 
energies. In other words, instead of the operator QoiE) which continuously depends on 
energy one obtains the discrete set of operators 

0g = ^/ go{E)dE, fc = l,...,iV, (50) 

Ek Jx)i, 

each of which is the averaged resolvent operator on the “on-shell” energy bin Tik where 
Dk is its energy width. 

4-2. Correspondence between “continuous” and “discrete” quantities 

Within the WP-formalism the discretization procedure involving three steps is intro¬ 
duced: 

(i) Division of continuous spectrum of the free Hamiltonian onto non-overlapping in¬ 
tervals and introduction of free WPs. 
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(ii) Projection of the scattering (and also bound-state ) wave functions and operators 
onto the above WP space. 

(iii) An additional energy-averaging procedure for energy dependent operators. 

It would be very useful to demonstrate how this discretization procedure works in prac¬ 
tical calculations by the example of solving the Lippmann-Schwinger equation for the 
transition operator t{E): 

t{E) = v{E)+vgo{E)t{E). (51) 

After application of the above three steps (i)-(iii) one gets a discrete set of operators 
at A G Dk (instead of continuous operator t{E)). The matrix elements of in the WP 
basis are related directly to off-shell elements of t-matrix: 

fife] / Q ^ \ 

g'GS, . (52) 

\ E J 

These operators satisfy to simple matrix equations 

1'= = 0 -f O0gi^ A € Sfc (53) 


where we denote by Gothic letters the WP projections of the respective operators. From 
eq. (ESI), one can get any on- and off-shell t-matrix elements whose energy and momentum 
dependencies are represented by histograms. It should be emphasized that f-matrix 
constructed in the WP representation satisfies exactly the unitarity relation [2^ . 

Let’s note that to find an elastic amplitude at energy E GT)k one has to solve eq. (ESI 
for one column of t-matrix only, i.e. i = 1,..., A^. 

Finally, the S'-matrix (and partial phase shift) can be found from the relation: 

^(A) « 1 - AeSfc, (54) 

Uk 


where Dk is the bin energy width. 

By similar derivation one can build the WP analogs for the Moller wave operators, 
total Hamiltonian resolvent, etc. (^ . 

The Table[T]shows the one-to-one correspondence between “discrete” and initial “con¬ 
tinuous” quantities. 

As a good illustration of such a fully discrete technique for finding the transition 
operator we present here the solution of a-a scattering problem where interaction includes 
both nuclear and Coulomb potentials 2l|. The basic s-, d- and g- partial a-a phase shifts 
found using the above WP-technique with additional energy averaging are displayed in 
Fig. El 


5. Multichannel scattering problem 

Very often many-body scattering problem in atomic, nuclear and molecular physics 
can be reduced to a multichannel scattering although the number of channels to be 
incorporated may be very large. Also a multichannel scattering problem arises in a 
simple two-body case when an interaction potential is not spherically symmetrical. 
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Table 1: Comparison between the continuous representation for the scattering theory basic objects (in 
the momentum space) and their discrete analogs in the WP subspace for E S T>k, q d'Zi and q' S Xtp 



Continuous 

Discrete WP 

1. The free resolvent 

go{E-q,q') 


2. The total resolvent 

9[E-,q, q') 

0M' ■ 

3. The t-matrix 

t{E-,q,q') 

tM'- 

4. The partial phase shift 

m 




Figure 6: s-, d- and g-partial phase shifts of o — a scattering found via the WP approach (solid curve) 
and from the direct solution of the Schroedinger equation by the Numerov method (dashed curve). 


The total multichannel Hamiltonian of the system can be written in a matrix form 
(in the channel indices v, v') 

hiyi/' — ~\~ Vi/i/' — 1 ,..., K ^ ( 55 ) 

where h^i, is the free-motion Hamiltonian in the channel v with simple continuous spec¬ 
trum oo), is the channel threshold, and are the coupling potentials. Further 
we will denote the operator matrices by the bold characters. The matrix of the free mul¬ 
tichannel Hamiltonian hg is diagonal and therefore the continuous spectrum of hg is the 
union of spectra of Hq^. Thus, the continuous spectrum of the multichannel Hamiltonian 
is degenerate (in contrast to a single-channel case), and the multiplicity of the degen¬ 
eracy r/ being equal to the number of open channels and hence depends upon energy. 
This means that at each energy there are rj scattering wave functions corresponding to 
different boundary conditions. 
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5.1. Free WP basis for a multichannel scattering 

Now one has to define free WP sets for each free motion channel Hamiltonian Hq^. 
We found that it is very convenient to make the bin partitions for continuous spectra in 
various channels so as the bins in different channels would be coinciding. In this case, the 
interval endpoints will be equal to the same energy values for all open channels. By this 
way one gets a degenerated discretized spectrum of the multichannel free Hamiltonian 
ho. Such a discretization method, as will be demonstrated below, makes it possible 
to determine multichannel t-matrix even without solving any scattering equations using 
properties of multichannel pseudostates [ 2 ^ . 






P.A' 






s: 


Q, n, 


Q, 


Q, 


hgj hg2 hgj 


Figure 7: The structure of the discretized spectra of the separate free motion Hamiltonians = 

1,2,3, (on the left) and the resulting degenerate dicretized spectrum of the total free Hamiltonian ho 
(on the right). The dashed lines represent bin boundaries while the solid lines inside bins correspond to 
free WP eigenvalues S*. 

So that, in the multichannel problem, one makes a partition of the total continnnm as 
usual; the region [H,, illmax] is divided into finite number of energy bins Si = [fi_i, 
Hereby an energy threshold should coincide with left end of some bin which we denote 
as * 1 /. Thus, for bins with indices i < iu there exist only u—1 open channels. We assume 
that channels are enumerated in the order of increasing their threshold energies and 
ii = 1. Figure [7] demonstrates how the dicretized three-channel spectrum is constructed. 

Using such partitions, the corresponding set of free WPs in each initial channel v. 

^ f dE\ipQ{E)), v = l,...,K, i = i^,...,N (56) 

V Hi JSi 

can be constructed from the exact wave functions Ii/jq^E)) of the free Hamiltonian 
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After the introduction of multichannel WPs it is straightforward to build a finite 
dimensional representation for all channel resolvents = [A + iO —completely 

analogously with a single-channel case. The channel resolvents enter the expansion for 
the multichannel free resolvent. So that, this operator has the following finite-dimensional 
approximation in the WP basis: 

K N 

0o(£^) = EEi<)5r(£^)«i- (57) 

y—1 i=ii, 

where complex eigenvalues are defined by the formulas similar to the single-channel ones. 

Further one can project the basic scattering operators onto such multichannel WP- 
basis similarly to the single-channel case treated above. In particular, the transition 
matrix elements for multichannel scattering problem is defined from the multichannel 
Lippmann-Schwinger equation: 

K 

{E) = E^ v,v' = (58) 

/i=i 

So that, a solution for the multichannel scattering problem can be found by solving the 
matrix WP analog of the above equation which is a direct generalization of the single¬ 
channel equation (15311 . 

However, there is another way by which one can solve the multichannel problem in 
the WP representation using properties of multichannel pseudostates. 


5.2. Eigenchannel representation and a multichannel resolvent 

To study the resolvent of the total multichannel Hamiltonian h, we shall employ the 
multichannel formalism in the so called Eigenchannel Representation (ER). The ER is the 
representation in which the multichannel S'-matrix takes a diagonal form with respect to 
the channel indices [s^ . More definitely, one can define at each energy E two orthogonal 
sets of the scattering functions for the Hamiltonian h: 

- wave functions {E))}^^j^ corresponding to the incoming waves in the channel 
iz (the so called experimental channel representation) and 

- wave-functions {|'*A^(£'))}^=i defined in the eigenchannel representation. 


The ER states differ from the experimental channel states by a rotation in the channel 
space with the rotation matrix dependent on the energy E @0. 

Just using the above ER formalism one can define multichannel scattering wave- 
packets and derive an analytical finite-dimensional representation for the total multi¬ 
channel resolvent (similarly to a one-channel case). 

It is convenient to construct multichannel scattering WP’s as integrals of the exact 
scattering wave functions for the total Hamiltonian h defined in the ER: 


\zk) 




wiE)\^-{E))dE, 


k = k ^,..., 


(59) 


where = [e^_ , e^] are new set of the total Hamiltonian bins whose parameters might 
be different for different >f. 
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Now we have to derive a spectral expansion of the total multi-channel resolvent using 
multi-channel scattering states defined in the ER. The spectral expansion of the resolvent 
in the “ordinary” set of the initial eigenfunctions of h can be written as sum of bound- 
state and continuum parts of the total resolvent, viz. g{E) = g^{E) + g'^{E), where 


Nt 

9^{E) = Y. 

rib = l 


\'4^nb){'4’nh I 

E - Erif, 


^ pOO 

= E / 


dE' 


\r{E')){riE')\ 

E + iO-E' 


(60) 


It is easy to show, that the following relation is valid: 

Y \r{E')){riE')\ = Y \riE'm-{E% (6i) 

U=1 >c=l 


which corresponds to a rotation between experimental channel and eigenchannel repre¬ 
sentations. 

Thus, we arrive at the spectral expansion of the continuum part of the total resolvent 
in the Hamiltonian eigenfunctions defined in ER: 


^ POO 

sbB) = E / 

>C=1 


dE' 


\^-iE')){i^-iE')\ 

E + iO-E' 


(62) 


where the thresholds in the ER coincide with thresholds of the initial experimental 
channels v. 

Applying the projection onto the multichannel basis (1591) one gets the following finite¬ 
dimensional representation for the total multichannel resolvent expressed via the multi¬ 
channel scattering WP basis 

K N 

= YY. \^X{E){z';:i (63) 

3pc— 1 k—k^ 


where complex eigenvalues g^{E) are defined as integrals over discretization bins, similar 
to the single-channel ones. 

Now the question arises, how to construct the states (15^ without solving the scatter¬ 
ing problem as we have done in a one-channel case. For this purpose we employ a new 
treatment of multichannel pseudostates. 


5.3. Multichannel pseudostates and a solution of the multichannel seattering problem 
without scattering equations 

The continuous spectrum of the total multichannel Hamiltonian h coincides with 
the spectrum of ho, so that at each energy there should exist g independent solutions 
corresponding to different boundary conditions for the scattering problem. Therefore, 
the usual L 2 discretization procedure becomes rather unclear here, because when the 
total Hamiltonian is represented by a respective matrix in some arbitrary finite T 2 -basis 
the required multiplicity of the spectrum disappears and one has only one pseudostate 
at every discrete energy value which can hardly be treated properly. So, it seems a one- 
to-one correspondence between the continuous spectrum (which is multiply degenerated) 
and its discretized analog is lost in the multichannel case. 
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Figure 8: The splitting of the eigenvalues in the degenerated discretized spectrum of the multichannel 
Hamiltonian ho caused by addition the interaction v. Qi/ are the channel thresholds while is the 
isolated eigenvalue corresponding to the single (in this example) bound state of the total Hamiltonian 
h. Sublevels inside each splitted set corresponding to the same eigenchannel number x are shown by the 
same lines (i.e. solid, dashed and dotted). 


However, the problem still can be solved if the discretized spectrum of the free multi¬ 
channel Hamiltonian is degenerate with necessary multiplicity as we constructed above. 
As was shown in ref. , the application of the perturbation v to the unperturbed opera¬ 
tor ho with degenerate discretized spectrum leads to a splitting of each ry-fold degenerate 
eigenvalue £* into the set of rj perturbed eigenvalues of the total Hamiltonian 

h as is schematically shown in Fig. [8] The very important property of the perturbed 
spectrum was shown in (on the basis of results from ref. [3l|), namely, these splitted 
sublevels form different intermittent sets in the energy domain between thresholds (i.e. 
the domain with fixed multiplicity rf) {e*^}i(LT ^ each of which corresponds to the xr-th 
branch of the continuum spectrum. So that the splitted levels are arranged in the same 
order (see Fig. 0): 

er+'>er, ^ = 

Moreover, this means that bands of splitted sublevels corresponding to different initial 
bins i are not overlapped with each other, which allows us to extract different k branches 
from multichannel pseudostate spectrum. 

Thus, one can select the new set of eigenvalues (and corresponding func¬ 

tions of pseudostates) for each value x which can be treated as the x-th branch of the 
discretized multichannel continuum and be confronted with the ER scattering states for 
respective branch of the continuous spectrum of the total Hamiltonian [^ . So, in this 
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way, we can establish a one-to-one correspondence between the discretized spectrum of 
the multichannel Hamiltonian matrix and the continuous spectrum of the initial multi¬ 
channel Hamiltonian. Due to such a classihcation of branches in the discretized spectrum, 
one can treat multichannel pseudostates on the same footing as the one-channel pseu¬ 
dostates. Now we can make the next step — consider multichannel pseudostates as 
approximations to multichannel WPs dehned in the ER (IKiB (we have done the same 
step in a one-channel case). 

In this way, we get the multichannel scattering WP which are expressed from the free 
WPs by a simple orthogonal transformation quite similarly to a one-channel case: 

rj N 

(64) 

u—l i=iu 


Now one can construct the WP approximation for the total multichannel resolvent via 
the relation (1551) and get a hnal solution of the multichannel scattering problem. 

As a good illustration of this approach we have calculated the coupled-channel S- 
matrix and the partial phase shifts and the mixing angle e for the Nijmegen I NN 
potential [s^ in coupled triplet ^Si — ^Di channels for the total spin S' = 1 and isospin 
1 = 0. The two-channel t-matrix has been obtained from the total resolvent using the 
well known formula 

t{E) = v-h vg(£')v, 

where finite-dimensional representation for the resolvent (16,SI) in the two-channel pseu¬ 
dostate basis is used. More definitely, the multichannel t-matrix for E G Di is defined 
by the relation: 


U.'{E) 


) , «|u|zb)(^bkK) 


V 




D, 


Dr{E - el) 


>C—1 k — kye 


{xl\v\z 
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k / ( rp\ 

9k \E) 


\v\x. 


/A 


(65) 


where | Zb) is a WP approximation for the bound-state wave function with energy and 
matrix elements of the interaction potential v between free and scattering WP states can 
be written using expansion ([M]) as follows: 

rj N 

{<\v\zk) = E E (66) 

^=1 3 = 3 , 1 . 

Thus, in the developed approach we can find the accurate multichannel t-matrix without 
solving any scattering equations at all. 

In Fig.[5]the partial phase shifts 6o{E), S 2 {E), and mixing parameter e{E) (i.e. those 
which enter the NN differential scattering cross sections - in the Stapp parametrization) 
are shown in very wide energy range 0 < Eiab < 800 MeV. They are found from just a 
single diagonalization of the respective two-channel Hamiltonian of the NN interaction 
using a two-channel free WP-basis. We compare these quantities in Figure to those ob¬ 
tained from direct numerical solutions of the two-channel Lippmann-Schwinger integral 
equation at many energies in the above energy range. 

This comparison shows that one can derive the accurate multichannel t-matrix or total 
resolvent for general multichannel case using a single diagonalization of the Hamiltonian 
matrix, i.e. without solving any scattering equations. It opens a new way in complicated 
coupled-channel calculations in atomic, molecular and nuclear physics. 
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Figure 9: The partial phase shifts 5o (S'): <52 (b) and the mixing parameter s (c) for the Nijmegen NN 
potential found from two-channel Hamiltonian diagonalization in the WP basis (solid curves) and from 
the direct numerical solution of the two-channel Lippmann-Schwinger equation (dashed curves). 


6. Stationary wave packets in a three-body problem 


Let’s consider now a general three-body scattering problem for particles 1,2 and 3, 
interacting via pairwise short-range potentials Va (a = 1,2,3). It is convenient to use 
here three Jacobi momentum sets (pa, Qa) corresponding to three channel Hamiltonians 
Ha {a = 1, 2, 3) which define the asymptotic states of the system. The symbol Pa denotes 
the relative momentum of b and c particles while qa is the momentum of the third particle 
relative to the c.m. of the pair {he}. In a general case of different particles the respective 
wave-packet bases should be constructed independently for each Jacobi set 22, 2^. Below 
we show how to built one of those. 

The channel Hamiltonian (e.g. Hi for a = 1) takes the form of a direct sum of two 
sub-Hamiltonians 

Hi = hi® hi, (67) 


where sub-Hamiltonian hi = hn ® vi includes the interaction vi in the {23}-subsystem 
and the sub-Hamiltonian hi corresponds to a free relative motion of the center of mass 
for particles 2 and 3 and the spectator particle 1. 


6.1. The lattice three-body basis 

The three-body wave-packet basis is constructed at first for the three-body free Hamil¬ 
tonian defined in the same Jacobi set: 

Ho = ho® hi. ( 68 ) 

So that the three-body free wave-packet states are built as direct products of the respec¬ 
tive two-body wave-packet states. To construct the three-body free WP basis functions, 
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we introduce partitions of the continua of two free sub-Hamiltonians, /iq and /iq, onto 
non-overlapping intervals {T>i = [£i-i,£i]}^i and {j)j = [£j-i,£j]}^^i respectively and 
introduce two-body free WPs as in eq. (051) : 



(69) 


where p and q are Jacobi momenta and Bi^Bj and f{p),f{q) are normalization factors 
and weight functions respectively. Here and below we denote functions and values corre¬ 
sponding to the variable q with additional bar mark to distinguish them from functions 
corresponding to the variable p. 

When constructing the three-body WP basis one should take into account spin and 
angular parts of the basis functions. Thus the three-body basis function can be written 
as: 



(70) 


where |a) is a spin-angular state for the {23} pair, \f3) is a spin-angular state for the 
third particle, and |r) is a set of three-body quantum numbers. The state (1701) is the 
WP analog of the exact free motion state in three-body continuum \p, q\a,f3 : P) for the 
three-body free Hamiltonian Hq [32|. 

The three-body free WP basis functions dZQl) are constant inside the cells of the 
momentum lattice built from two one-dimensional cells and momen¬ 

tum space. We refer to the free WP basis as a lattiee basis and denote the respective 
two-dimensional bins (i.e. the lattice cells) by 'Dij = TiiiSi 

Such multi-dimensional wave-packet basis is an “eigenbasis” for the free Hamiltonian 
Hq and has the same properties with respect to this operator as two-body free WPs 
with respect to the two-body free Hamiltonian Hq. In particular, every operator which 
functionally depends on Hq has an explicit finite-dimensional representation in the lattice 
basis discussed (^ . So that, it is straightforward to obtain a diagonal representation 
for the three-body free resolvent Go and to use it further for solving few-body scattering 
equations. This is a way to develop the wave-packet matrix scheme which is a finite¬ 
dimensional form of the initial scattering equations. 

However, the above discrete representation has a remarkable advantage in comparison 
with the continuous one: one can construct and employ a few-body wave-packet basis 
directly for the three-body channel Hamiltonian Hi. This allows to simplify a solution 
of scattering problems drastically. 

6.2. The WP basis for the channel three-hody Hamiltonian 

Let’s introduce a WP basis for two-body sub-Hamiltonian hi. Assume as above that 
there are N\y bound states in the {23} subsystem with corresponding bound-state wave 
functions {\zf^)}^2i and eigenenergies {e“*})))bi. One defines the partition {^k}k=Ni,+i 
of the continuous spectrum of hi and constructs the set of scattering wave packets \z^) 
from the respective exact scattering wave functions |'0p) according to eq. (H51) . The com¬ 
plete WP basis {\z'^)}k2i foi' the hi sub-Hamiltonian includes bound state functions and 
scattering wave packets and its functions may depend on possible spin-angular quantum 
numbers. 
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Now one can build the three-body wave-packets (3WP) for the channel Hamiltonian 
Hi just as products of two types of WP states for hi and /ij sub-Hamiltonians whose spin- 
angular parts are combined to the respecive three-body states with quantum numbers 


P: 


\Zlf) = \z^,x^,a,l3:r), = j = l,...,N. (71) 

The properties of the 3WP constructed in this way are the same as properties of the 
two-body wave packets viz. they form an orthonormal set and any operator functionally 
dependent on the channel Hamiltonian Hi has a diagonal matrix representation in the 
subspace spanned on this basis. It allows us to construct a finite-dimensional approxima¬ 
tion for the three-body channel resolvent Gi{E) = [i5 -|- iO — Hi]~^ in a very convenient 
analytical form. 

Indeed, the exact three-body channel resolvent is a convolution of the two-body sub¬ 
resolvents gi{E) = [E + i0 — hi\~^ and gl{E) = [E -\-iO — hl]~^ 22, 23, 13, 


Gi{E) 


1 

27ri 


de5i(e)ffo(£’ 


e). 


(72) 


Using spectral expansions for these two-body resolvents and making the integration, one 
gets an explicit expression for the exact channel resolvent Gi as a sum of two terms 
Gi{E) = Gf'^^E) + Gi'", where the bound-continuum (BC) part takes the form [2^ : 


G 


BC 


Nb 

(^) = E U 

r.a,^ra=l' 


dq 




E + iO-e^* - 


SL 

2M 


(73) 


where M is the reduced mass in the {23} -I-1 channel. While the continuum-continuum 
(CC) part of Gi takes the form: 



dg 




E 


. iO - Pi - 

2M 


r| 


(74) 


where g is the reduced mass in the {23} subsystem. 

Projecting further the exact channel resolvent onto the three-body channel 3WP basis 
defined in eq. CD, one can find analytical formulas for the matrix elements of the Gi 
operator. The respective matrix is diagonal in all wave-packet indices: 

= E E \Zlf)Glf{E){Zlf\. (75) 

r,a,^ kj 


The matrix elements G^“^ are defined as integrals over the respective momentum bins 
for the BC part: 


_ * 

Cri,, — 


\f{q)\^dq 


kj 


B, 


U -b iO - e“* - 


■ 

2M 


and for the CC-one: 


qToiP _ ^ 



f |w(p)P|/(g)pdpdg 
4 U + iO-g-^ 


k=l,...,Nh, j = l,...,N 

, fc = TVb-b I,... ,iV“, j = 


CS*) 


I,...,7V. 
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These elements depend on the spectral partition parameters (i.e. and 'Dj bin end¬ 
points) and total energy E only. They do not depend explicitly on the interaction 
potential vi. When solving the scattering equations in the finite-dimensional WP basis 
the corresponding solution converges with increasing the basis dimension, the final result 
turns out to be independent of the particular spectral partition parameters. Explicit 
formulas for the resolvent matrix elements (175b ') and (I75b 'l are given in Appendix A. 

The explicit analytical representation (123 for the channel three-body resolvent is a 
basic feature for the wave-packet approach since it allows us to simplify solution of the 
general three-body scattering problem drastically. In particular, this representation has 
been used to solve the finite-dimensional analog for the Faddeev equations 


27, 


[. This 

representation has been employed also to solve some particular three-body scattering 
problems using the three-body Lippmann-Schwinger equation 2^, 2^ ^ for composite 
projectile scattering off nuclear target (see below). 

To use the states dZID practically, one can approximate them with the pseudostates 
of the sub-Hamiltonian hi in some L 2 basis. As has been shown earlier in Section 3, the 
free WP basis is very appropriate to approximate scattering states because the respective 
functions have a very long-range behavior in configuration space which makes it possible 
to describe properly an asymptotical behavior of Faddeev wave function components at 
large distances. So one can calculate the eigenstates (the bound and pseudostates) of 
the sub-Hamiltonian hi matrix in the two-body WP-basis {|a;i)}i^i via a diagonalization 
procedure. As a result one gets the eigenstates of the hi sub-Hamiltonian expanded in 
the free WP basis (for each quantum number a) similar to eq. (H!?)) . 

Hence, starting from the free WP bases for each two-body sub-Hamiltonian one gets 
a set of three-body basis states both for free and channel Hamiltonians, Hq and Hi 
respectively, which are related to each other by a simple matrix rotation. Using projection 
of operators and wave functions onto these 3WP bases it is possible to solve three-body 
scattering problems via matrix formalism. 

It should be mentioned that in a general few-body case, the free and the channel 
wave-packet bases can be constructed by a very similar procedure using a continuum 
discretization for every subsystem. In the bases constructed in this way, it is straight¬ 
forward to find explicit finite-dimensional representations for channel and free resolvents 
and furthermore to solve the resulting scattering equations in a very convenient matrix 
form. 


7. Composite particle scattering oflf a nuclear target 

Let us consider a few-body problem of an elastic composite particle scattering off 
a target when rearrangement channels are neglected while intermediate excitations into 
continuum are taken into account. 

7.1. Solution of the Lippmann-Schwinger equation 

The total Hamiltonian for the composite projectile plus target system can be written 
as follows: 

H = hint + he + 14xt, (76) 

where hint = hg + [rjj ) is the internal sub-Hamiltonian for the projectile com- 

i<j 

posed from several fragments i = 1,...,A", while hc{R) is the sub-Hamiltonian of 
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projectile asymptotic motion including the center of mass Coulomb interaction, and 
14xt = are the fragment-target external interactions where are separate frag- 

i 

ment positions and are their relative distances. 

The optical potentials Vi are usually assumed to be complex and energy-dependent. 
However, the smooth energy dependence of the input optical potential is of no significance 
for our present purpose and, thereby we omit it here. Also we assume that there are 
bound states in hint sub-Hamiltonian {\zn)}nLi with energies 

The channel Hamiltonian is a direct sum of two sub-Hamiltonians: 

Hch — hint ©he. (77) 

In the case to be discussed, when rearrangement channels are neglected, the transition 
operator T which describes the elastic scattering of the projectile by the target nucleus 
as well as the projectile breakup is found from a single Lippmann-Schwinger equation 
(LSE) [H: 

T{E) = Hext + HextGch(A)T(£;), (78) 

where Gch = [E + iQ — Hch]~^ is the channel resolvent which determines the asymptotic 
states in the elastic channel. 

Then, wave-packet basis states corresponding to a discretization of the channel Hamil¬ 
tonian continuum are built as products of hint and he sub-Hamiltonian WP basis states 
as defined in eq. (EB). The only difference in that we have to replace free WPs over the 
center of mass variable with Coulomb WPs for the sub-Hamiltonian he- Thus our WP 
basis in this case takes the form: 

|Z4'^) = |z[,xf;/,A:T), (79) 

where I and A are subsystem orbital momenta while A is a total orbital momenta of the 
system which assumed to be conserved. 

Thus, after a WP-projection of wave functions and scattering operators one gets 
respective vectors and matrices in the channel WP space. Such a wave-packet represen¬ 
tation is eigen for the projected channel Hamiltonian and the basis states from eq. EH) 
correspond to exact asymptotic states of the system. In particular, the initial state 
wave function of the system which defines the asymptotic free motion of the projectile 
corresponds to a single WP-state. For example, if one study an elastic scattering of 
the projectile in its ground state (assume that 1 = 0), the initial state wave function 
|d>o (A)) = |z°, ijj^^E — e*); L) is represented by the state of the channel WP basis 

where A — ej € Ttjg. 

Eventually, after the wave-packet projecting of the scattering operators, all terms in 
the LSE (17811 are reduced to a matrix form and thus the T-matrix can be found from a 
single matrix equation: 

T = V + VGchT, (80) 

where V is the external interaction matrix in the three-body wave-packet basis and Gch 
is a diagonal matrix of the channel resolvent E3 taken at the total energy E. It should 
be stressed that to find on-shell and half-shell T-matrix elements, it is sufficient to solve 
the matrix equation (18011 only for one column = T„_„c, n = , N ■ N (we use here 
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Figure 10: The comparison of the deuteron off ^®Ni-target elastic cross sections (in ratio to the Rutherford 
cross section crji{9)) obtained in the WPCD approach with an inclusion of intermediate deuteron breakup 
channels (full curve), with the folding model, i.e. without taking into account intermediate breakup 
channels, (dashed curve) and the experimental data (circles) [33| . 


the multi-index n = (A:, j), so that the label no = (1, jo) corresponds to the initial state). 
The respective solution can be found from the equation 


(I - YGch)t = V. 


(81) 


Here the notation is: Vn = n = 1,..., N■ N, and I is the unit matrix. Finally, the 

S-matrix elements are interrelated to the vector of t-matrix elements by the relationship 


^jo 


(82) 


As a particular example for such a scattering problem we consider the well known test 
case: the deuteron elastic scattering off ®®Ni target — see the respective results in Fig. [TUI 
The details of a calculation and potential parameters used can be found in ref. 2^ . Also 
we have shown [2^ that our discrete WP-solution for such a coupled-channel problem 
is in a very good agreement with that derived from the traditional CDCC approach in 
which one has to solve a large system of coupled differential equations in every partial 
wave. 


1.2. Construction of an effective projectile-target interaction 

In many problems of quantum physics it is necessary to know an operator of effective 
interaction between a composite projectile and stable target or vice versa, which takes 
into account the virtual excitation of the incident particle (or target), i.e. its dynamic 
polarization in the scattering process. Such an operator of the effective interaction is 
replaced usually by some phenomenological optical potential between composite projec¬ 
tile and stable target (or vise versa) because its microscopic evaluation is often a very 
cumbersome problem. So the convenient approach for a theoretical calculation of the 
effective operator is a very good object for the theory. 
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The wave-packet discretization method allows us to construct such operators explic¬ 
itly, using the finite-dimensional representation of the channel resolvent obtained above. 
In this solution it is convenient to apply the Feshbach projection operator formalism . 

For the reader’s convenience, we briefly recall a derivation of the effective operator of 
interaction within such projection approach for the case of composite particle scattering 
off a nucleus discussed in the previous subsection. Consider the Schroedinger equation 

H\-^) = E\^) 

for the wave function of the system defined by the Hamiltonian (1751) . 

Let, l^i) be the projectile ground state with energy e^. Using the above projection 
operator technique, one defines two projection operators: 

F=\zi){zil Q = l-F, QF = 0. (83) 

Here F is the projector onto the elastic channel, Q is the orthogonal projector onto all 
the inelastic channels. As a result the Schroedinger equation for the total wave function 
is splitted into two coupled equations for components F\'if) and Q|'I') : 

{FHF - F)F\'ii) = -FHQl'S), 

{QHQ-F)Q\'if) = -QHF\-^). (84) 

Further, we introduce the resolvent of the Hamiltonian projected onto Q-subspace: 
gI^\f) = [{F -I- iO) — QHQ]~^ . Then, substituting the second equation in (IMll to 
the first one we obtain the well known equation for the elastic scattering function Ajd'): 

{FHF + FUes{E)F - F)F\'^) =0, (85) 

where 

FUes{E)F = FHQG^^^QHF ( 86 ) 

is an effective nonlocal and energy-dependent operator of interaction between a composite 
particle and the target with taking into account all the inelastic channels. Since the 
operators F and h-mt 0 he commute, and the orthogonality condition FQ = 0 is valid, 
the equation (1851) takes the form of a conventional two-particle Schroedinger equation for 
projection |xi) = (zild*) of the total wave function onto the elastic channel: 

[he + Vfoid + C/eff -{E- el)]\xi) = 0. (87) 

Here we have introduced a folding potential Ufoid describing the interaction of the pro¬ 
jectile center of mass with the target nucleus when the inelastic channels corresponding 
to the excitation or breakup of the incident composite particle are fully neglected: 

VfoM{R) = {zi\Ve^t\zi). (88) 

Here the integration in the the matrix element is carried out only over the internal 
variables. All inelastic (and breakup) effects are included into the non-local Feshbach 
potential: 


Uee{E) = (zilUextQGQ(A)gUextl^l). 
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(89) 


Despite the compact notation, the computation of this operator is in principle even more 
difficult problem than the solution of the original scattering problem. However, when 
using the wave-packet basis (17^ consisted from eigenfunctions of the channel Hamiltonian 
the practical construction of operator Ueg{E) becomes clear and straightforward. 

Indeed, using the wave-packet states for the channel Hamiltonian, one finds the finite¬ 
dimensional diagonal representation for the channel resolvent Gch{E) obtained via for¬ 
mulas from eq. ( 123 . 

In the WP-approach the Q-subspace is just orthogonal part of the WP subspace, 
so that the finite-dimensional representation for the projector Q can be written as (we 
omitted possible spin-angular indices): 

£l = ^\zk,xf){zk,xf\, (90) 

k,j 


where upper prime symbol means that the sum does not include the ground state of hint- 
Thus the matrix representation for the projected Gq operator is straightforward: 

^QiE) = mliE)]-^ - ( 91 ) 

where all operators are meant as being projected onto the Q-subspace. 

After direct evaluation of the matrix for GQ-operator from eq. (lOTl) . the explicit 
formula for the effective potential in the wave-packet representation can be written in 
the form: ^ ^ 

C/eff(A, R, R')=J2J2 Bkjm<GQiE)]k,,k'j'Bl,^,iR'), (92) 

kj k'j' 

where the form-factors are defined as the integrals: 

Bfc,(A) = (zi|Hext|^fe,xf), (93) 

and an integration is done over internal variables only. In practical calculations, all the 
necessary spin-angular parts of wave functions should be taken into account. 

Thus, the WP approach gives a direct and convenient way to calculate an effective 
optical potential for an interaction of a composite particle with a stable target, also this 
formalism can be generalized for constructing an effective interaction between two col¬ 
liding composite particles which may be excited or disintegrate in the scattering process. 

In Fig. [11] the real and imaginary parts of the effective nonlocal optical potential for 
the deuteron and the ®®Ni nucleus interaction at incident deuteron energy Ed = 80 MeV 
calculated via the WP-approach are displayed for the total orbital angular momentum 
L = 0. 

Besides the direct calculation of the complicated nonlocal interaction operator, the 
approach described above is very convenient in those cases where inelastic channels play 
a role of a correction to the elastic scattering when the main contribution comes from 
the folding potential. In such a case, one can employ the “inner” few-body WP basis of 
a rather small dimension for the calculation of the effective potential and the “external” 
two-body WP basis of a large dimension for the subsequent solution for the resulting 
two-body scattering problem with the total interaction including the folding potential 
and the above effective potential (2^ . 
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Figure 11: The real (left) and imaginary (right) parts of the effective deutron-®®Ni interaction optical 
potential calculated at -Ec.m. =80 MeV for the total angular momentum L = 0. 


Just such a situation arises very often when the energy of collision of a composite 
particle with a stable target increases. The relative weight of the breakup channels 
decreases (after some characteristic energy) and thus the contribution of the effective 
potential taking into account the projectile breakup is reduced as in our above example. 

As an illustration for such an approach we present in FiglT^] the differential cross 
section for elastic deuteron scattering off the nucleus ®®Ni at energy Ed =80 MeV calcu¬ 
lated using the above Feshbach potential. The result of such optical model calculation 
nearly coincide with the solution of a direct three-body problem found with the CDCC 
approach in ref. [1^. The dimensions of the Coulomb WP-basis |xf) used in Q projector 
is in five times less than the dimension of the same Coulomb WP-basis used in “external” 
projector F. 



„ (deg) 


Figure 12: The differential cross section for the elastic deutron-®®Ni scattering at Bc.m. =80 MeV 
calculated in the CDCC approach (dash-dotted curve) ISHl and with an effective potential of the Feshbach 
type constructed in the WP approach (solid curve). 
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Let us summarize, with the described technique we are able rather easily to calculate 
the theoretical non-local optical potentials for a scattering of composite projectiles off 
stable targets (or vice versa), or even for a collision of two composite particles. 


8. Formulation and solution of the Faddeev equations in the wave-packet 
representation 

The most rigorous and correct formulation for the few-body scattering problem can 
be attained using the Faddeev or Faddeev-Yakubovsky equations. We will demonstrate 
below that the wave-packet approach being applied to solving such equations gives the 
great advantages. 

8.1. Formulation of Faddeev equations in momentum space and peculiar properties of 
their solution 

Here we consider a solution of the Faddeev equations for a scattering of three identical 
particles 1, 2 and 3 with mass m (nucleons). In this case the elastic scattering observables 
can be found from a single Faddeev equation (FE) for the transition operator U (the so- 
called AGS equation), e.g. in the following form [32|: 

U = + PtGoU, (94) 

where t is two-body off-shell t-matrix in three-body space, Gq = {E+iO—Ho)~^ is the free 
three-body resolvent and P is the permutation operator which changes the momentum 
variables from one Jacobi set to another one. For the case of three identical particles the 
operator P is defined as sum of two cyclic permutations of particles: 


P = P 12 P 23 + Pl3P23- (95) 

It should be emphasized that a similar permutation operator is included in kernels of the 
Faddeev equations in the case of various particles. 

After the spin-angular expansion, the operator equation (IMl) for each value of the 
total angular momentum and parity is reduced to a system of two-dimensional integral 
equations in momentum space (or to coupled two-dimensional integro-differential equa¬ 
tions in the conhguration space). Although Faddeev [s^ proved that the kernels of these 
equations belong to the Fredholm type (i.e. the inhomogeneous equation has a unique 
solution), a practical solution of coupled two-dimensional integral equations is a very 
complicated and time-consuming task due to complicated singular structure of integral 
kernels and also the large number of coupled spin-orbital channels. 

One of problems is that kernels contain two-body off-shell t-matrices t{q,q',E) for 
different values of the total orbital angular momentum and spin of the pair interacting 
particles, which should be calculated many times at different energies. 

The second problem here is that the Faddeev-type kernel at real energy contains sin¬ 
gularities of two types: two-particle cuts, corresponding to the presence of bound states 
in the subsystems, and the three-body logarithmic singularity (at the breakup threshold). 
The two-body singularities are easily eliminated by the technique of residues, while for a 
regularization of the three-body cut a number of special techniques have been suggested 
in the previous years. As a result, the whole solution procedure becomes rather compli¬ 
cated. Among such specific techniques the following are employed most often: 
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- choice of a special quadrature grid points for a momenta q' defined in one Jacobi set 
depending on a momentum q defined in another Jacobi set together with the spline in¬ 
terpolation of an unknown function at iteration method of solving; 

- solution of the equations at the complex energy plane followed by an analytic contin¬ 
uation to the real axis; 

- shift for an integration contour from the real axis to the complex momentum plane. 

However, the main characteristic feature of the Faddeev kernels is the presence of the 
particle permutation operator P. The kernel of this operator P{p, q;p', q') as a function 
of the momenta contains J-function and two 0-functions and this results in variable 
integration limits for integrals in the Faddeev kernels (j^ . Therefore, when replacing the 
integrals with respective quadrature sums in numerical procedure it is necessary to apply 
an interpolation of the unknown solution depending on several variables at each step of 
an iterative procedure. As a result of such elaborated multi-dimensional interpolation on 
each iteration step, the interpolation procedure takes most of the computational time of 
the whole computation and requires usage of powerful supercomputers Q)- 

On contrast to this, the WP method described here allows to circumvent completely 
the above difficulties in solving the Faddeev (and Faddeev-Yakubovsky) equations. First 
of all, instead of eq. (1941). one uses the equivalent form of the Faddeev equation for the 
transition operator U |29l. 1^ : 

U = Pvi + PviGiU. (96) 

In eq. (IMl) Gi is the three-body channel resolvent introduced in Section 6. Due to the 
identity tGo = viGi the kernels of the equations (1941) and (IMl) are the same and their 
solutions U and U coincide on-shell and half-shell. 

This form of the equation is especially useful in the WP representation with basis 
functions corresponding to the channel Hamiltonian iJi, because in this representation 
the resolvent Gi has a simple analytical form and explicitly depends only on the partitions 
of the continuous spectrum of the subsystem. Thereby the need for multiple calculations 
and interpolation of the off-shell f-matrix is eliminated. 

Further, all the singularities of the Faddeev kernel in the form (1961) are concentrated 
in the channel resolvent Gi and they are smoothed when projecting on the WP basis (av¬ 
eraged by the integration over the energy bins). Therefore the resulting matrix equation 
can be solved directly for real energies. 

Finally, the use of matrix of permutation in the WP basis (as well as in any other fixed 
basis) completely eliminates the need for very numerous interpolations of the required 
solution at each iteration. 

So, all these innovations taken together lead to great simplification in practical solving 
the Faddeev-type few-body scattering equations. 

8.2. Three-body wave-packet basis for 5N system 

To illustrate this novel technique we consider the realistic Nd scattering problem, 
where one needs to particularize the three-body wave-packet basis for the case of three- 
nucleon system with tensor NN interactions. We use the following quantum numbers 
for the subsystems defined in Section [B] and in the whole three-body system according to 
the (jj)-coupling scheme: 

O: S, J23}; P : 
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r = {j,7r,r}. 


(97) 


where I, s and J 23 are the NN subsystem quantum numbers: I is an orbital momentum, 
s is a spin and j 23 = 1 + s is a total angular momentum of the subsystem (the interaction 
potential depends, in general, on the value of j 23 )- The other quantum numbers are the 
following: A is an orbital momentum and I = A + cr is a total momentum of the third 
nucleon, where cr = ^ is its spin. Finally, J = j 23 + 1 is a total angular momentum of the 
three-body system, T is a total isospin and tt is parity, all of them are being conserved. 
Let’s also note that the pair isospin t can be defined by values of I and s, because the 
sum I + s + t must be odd. 

The two-body free WP states (1551) should be defined for each partial wave I and A and 
further they are multiplied by appropriate spin-angular states according to dnH). Thus 
the free three-body basis function ([701 can be written in the detailed form: 

K“^) = k,s-;«,/3:r), (98) 

The three-body WP states corresponding to the channel Hamiltonian Hi are defined 
here as has been detailed in Section [0] However now one has to take into account all 
possible spin-angular couplings in {23} subsystem induced by tensor couplings in vi 
interaction, so the two-body scattering WPs of hi sub-Hamiltonian must be defined in 
the eigenchannel representation as in eq. dSO]). Then, three-body channel WP states have 
the following form: 

\Zlf^) = \z^,x^;a,P:r), (99) 

where the spin-angular quantum numbers a in {23} subsystem corresponds to the eigen¬ 
channel representation and differs from those, i.e. a, for the free WP basis. In particular, 
instead of an angular momentum value I this index contains eigenchannel value so 
a = {>c, 5 ,^ 23 } while a = {lsj 23 }. 

These states are the WP analogs for the three-body scattering states , q]a,j3 : P) 
of the channel Hamiltonian Hi, where \ipp) is a scattering wave function of hi sub- 
Hamiltonian defined in the ER. 

In Section 5, it was shown that the pseudostates obtained via the diagonalization 
sub-Hamiltonian hi in a free WP basis are good approximations for the exact WPs \z^) 
constructed from the scattering wave functions for hi. The exact WPs are related to the 
free WPs by a simple orthogonal transformation: 

N 

(100) 

i=l I 

where spin-angular parts of wave functions are taken into account as well and the multi¬ 
index a = {>c, 5 ,^ 23 } related to the ER is used. 

Now we have in our disposal the WP basis for the channel three-body Hamiltonian 
and can apply explicit formulas (1751) for the channel resolvent Gi which enters eq. (1551) . 

8.3. Matrix analog for the Faddeev equation in the WP basis 

In our approach, all the operators in eq. (1961) are projected onto 3WP basis corre¬ 
sponding to the channel Hamiltonian Hi. In other words, every operator, e.g. U, is 
replaced with its finite-dimensional WP representation: 

E E (101) 

kj,k',j' 
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As a result, one gets the matrix analog for the Faddeev equation (IMl) (for each value of 

r) 

U = PVi + PViGiU. (102) 

Here Vi and Gi are the matrices of the pairwise interaction and the channel resolvent 
respectively, the matrix elements of which can be found in an explicit form. 

Thus, to find the elastic scattering amplitude it is required: 1) to calculate matrix 
elements of P, Vi, Gi matrices and 2) to solve the system of algebraic equations (11021) . 

The matrix Vi of the potential vi is diagonal in the indices of the wave-packet 
basis for the free sub-Hamiltonian hg and has the block form: 

= Sp/S'Sjj'Sss' {Zk\vi\z^')333 ( 103 ) 

for a = {>!:, s, _) 23 }. Here the subindex j 23 in the potential vi matrix element means that 
it depends on the two-body total angular momentum value. The matrix elements (11031) 
do not depend on index j and can be written with the usage of the rotation matrix O 
dehned in eq. (llOOl) as: 

{zm4)333 = 


where {xl\vi\xl ,)are the potential matrix elements in the free WP basis. 

The matrix of the operator P in the free three-body packet basis corresponds to the 
overlap between basis functions defined in different Jacobi sets: 


rpOia/J.a'/?' 




{Xf4{l)\x4^'{2)) 


{Xf4 {1)1X14'(3)), 


(104) 


where the argument 1 (or 2 and 3) in the basis functions denotes a corresponding Jacobi 
set. Such matrix elements can be calculated by integration over the basis functions in 
momentum space: 


rpOia^.a'^' 



dp' dq' X 




did-i/ dj dj 


(105) 


where the prime at the lattice cell indicates that the cell belongs to the other 
(rotated) Jacobi set while the 4/3 a'P'iP^djP'd') is the kernel of particle permutation 
operator in a momentum space. This kernel, as is mentioned above, is proportional 
to the product of a Dirac delta and Heaviside theta functions. However, due to the 
integration over momentum bins in eq. (fT05l) these singularities are averaged over the 
momentum lattice cells and, as a result, the elements of the permutation operator matrix 
in the WP basis are finite and non-singular. 

The detailed technique for the calculation of the matrix element in eq. p05l) is pre¬ 
sented in ref. ^|. 

The permutation operator matrix P in the channel 3WP basis is expressed through 
the overlap matrix P° for the lattice basis function of eq. (11051) with help of the rotation 
matrices O if one uses the pseudostate approximation (HOOD for the scattering WPs |z|): 


{zl4\P\z4'^') 


ii' LI' 


r)>' 

ki ^ k'i' 




(106) 


So, we have a relatively simple formulas and numerical algorithms to determine all 
the quantities entering the kernel of the matrix Faddeev eauation p02D . 
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8.4- Determination of observables in the WP approach 



Figure 13: The differential cross section for the elastic nd scattering calculated via the WPCD approach 
with the Nijmegen I NN potential for incident neutron energies En = 3 MeV (solid curve) and En = 
9 MeV (dashed curve) in comparison with the experimental data (circles) [3^ . 


The elastic on-shell amplitude in the wave-packet representation is calculated as a 
diagonal (on-shell) matrix element of U-matrix [29| : 




2m 

3go 


(107) 




where m is the nucleon mass, go is the initial two-body momentum and 

; ao, /3 : r) is the 3WP basis state corresponding to the initial scatter- 


30 

ing state. Here is the bound state of the NN pair in the initial state (the deuteron, 
in our case), the index jo denotes the bin Ttjg including the on-shell momentum go and 
djg is a momentum width of this bin. 

In Fig. [T3] the differential cross sections for nd elastic scattering found with the Ni¬ 
jmegen INN potential in the WP approach are represented for different incident 
neutron energies in comparison with experimental data 37| (circles). It is evident from 
the Figure that agreement with the data is very well. 

In Fig. [U] the comparison is given for the neutron vector analyzing powers Ay for the 
elastic nd scattering at 35 MeV. Here the WP basis with dimension N x N = 100 x 100 
has been used and the partial waves with the total angular momentum up to J < 17/2 
have been taken into account. 
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Figure 14: The elastic nd scattering differential cross section (a) and the neutron vector analyzing 
power Ay (b) at 35 MeV obtained within the WP approach (solid line). The experimental sets are the 
following: pd data at 35 MeV (full circles), nd data at 36 MeV [3^1 (empty circles) and nd data 35 
MeV (triangles) correspondingly . 


8.5. Treatment of the three-body breakup in the WP-approaeh 

Now let’s pass to a treatment for the three-body breakup using the WP approach. 
One can show that the breakup amplitude may be defined as a matrix element of the 
same transition operator U satisfying eq. dMl), not only as its diagonal element on initial 
states, but also for the transition to the continuum state of the channel Hamiltonian Hi: 

( 108 ) 

pqqo 

where is the scattering function for the Hamiltonian hi corresponding to the 

outgoing boundary condition, p,q are the final momenta of the subsystem {23} and the 
third nucleon correspondingly, qq is the initial momentum of the third nucleon which are 

interrelated to each other by the energy conservation e} + ^ 

Thus, in the WP-approach, the breakup amplitudes can be defined quite similarly to a 
matrix element for the elastic scattering transition operator U with evident replacement 
of the NN bound-state wavefunction with the exact scattering functions for the NN 
sub-Hamiltonian (or the corresponding WP state). 

As an illustration, consider the case of semi-realistic s-wave NN interactions MT HI. 
In this calculation one takes I = X = 0 and one has only NN spin quantum number 
s to distinguish different spin-angular channels: ag = s = 1 for the initial channel, 
ex = s = 0,1 for the final channel and the three body quantum number P is defined 
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by the total spin value S. We show here the results for the hyperspherical breakup 
amplitude which defines an asymptotic behavior of a Faddeev component of the wave 
function and is related to the breakup amplitude (11081) as follows: 


A^{e) = 9 = arctan(^) 


3^3 


2p 


(109) 


where 9 is the hyperangle in momentum space. In the WP approach the breakup ampli¬ 
tude T is defined by non-diagonal elements of the transition operator il (see details in 
ref. [H): 




T^^P.q) 






tpSs 

JS(pI) lio-fcj 

PlQjqo ’ 

pe A*, 


( 110 ) 


where |^y|,) is the WP basis state corresponding to the initial state: index 1 denotes 
the bound state of the NN pair (deuteron) and index jo denotes the “on-shell” q-bin 
'Sjg, while \Z'^f) defines the state of the channel Hamltonian Hi “excitation” to which 
corresponds to a breakup process. Here 6(pi) is the s-wave phase shift corresponding to 
the NN pair interaction at energy ej®, pi = and q* = ^[qj-i + qj] are momenta 

corresponding to and Sj bins respectively and the dl is the momentum width of the 
Al bin. 

A comparison for the hyperspherical breakup amplitudes defined in the WP approach 
with results of the benchmark solution of the Faddeev equation is displayed in Fig.lTKl 


8.6. Some features of our algorithm for solution of Faddeev equations in the WP ap¬ 
proach. Realization through the GPU parallel computations 

Using the WPCD approach, we have reduced the practical solution of the Faddeev 
equation to solving the system of algebraic equations (11021) . Thus, the WPCD algorithm 
for solving system consists of the following main steps: 

1. Construction of free WP bases {|a;i), 1%)}, calculation of the potential matrix 
{x°j\vi\x°j ), diagonalization of two-particle sub-Hamiltonian matrices for each value 
of I and s and finding the energies of pseudostates and their functions in the WP 
representation, i.e. matrices 

2. Calculation of the permutation matrix P° in the lattice WP basis (11051) . 

3. Calculation of the channel resolvent matrix Gi{E) in the WP basis corresponding 
to the channel Hamiltonian Hi (the diagonal matrix). 

4. Solution of the system of algebraic equations (11021) and the determination of the 
elastic and breakup amplitudes. 

Point 1 provides two-particle input for a solution of the three-body problem. In our 
WP approach, the input is obtained by a single (for each value of a) diagonalization of 
the two-particle Hamiltonian. The result of these diagonalizations can be used to solve 
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Figure 15: A comparison of the hyperspherical breakup amplitudes A for the spin-quartet (a) and spin- 
doublet (b) channels obtained within WP technique (solid curves) and in the conventional Faddeev 
calculations |41|| (full circles). 


the scattering problem at different energies E. Whereas in the standard approach to 
solve the Faddeev equations one has to calculate two-particle off-shell t-matrices many 
times for each value of E. 

Point 2 is the key to simplify the whole solution of the Faddeev-type equations. The 
use of a finite-dimensional approximation for the permutation operator, i.e. replacing it 
with a fixed matrix (in any basis), avoids the need for numerous and time-consuming 
multi-dimensional interpolations of a current solution during the iteration process. In 
the standard approaches, these interpolations take the most of the computing time. Al¬ 
though in the calculation of matrix elements of permutation operator in our approach 
also meets some difficulties — multiple integrals with variable limits, however these 
matrix elements are computed with simple functions. Nevertheless, the calculation of 
the permutation matrix takes the major part of computing time in our algorithm 
in sequential execution on the CPU. But it should be stressed that the matrix is 
independent of energy and therefore being calculated at once, it can be used to solve 
scattering problem at so many energies as we wish. Note also that due to energy con¬ 
servation {p\/m + iq\/{Am) = p\/711 + where subscripts 1 and 2 indicate the 

different Jacobi sets, the permutations matrix is very sparse - only about 1 % its 
matrix elements are nonzero. 

Point 3 is based on the main feature of the WP representation. It is just the ad¬ 
vantage of the WP basis that the resolvent matrix Gi{E) for the channel Hamiltonian 
Hi is diagonal, does not depend explicitly on interaction potentials and its elements are 
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determined only by the partitions of two-dimensional bins and the total energy E using 
explicit formulas. 

The only difficulty in the practical solution of the matrix equation (I102D is its high 
dimensionality — that is the price for the elimination of the basic difficulties of a standard 
approach to solving the Faddeev equation. Solution of an algebraic system of such high 
dimensionality is an untrivial problem even for a supercomputer. But in our approach 
we do not need in solution for the whole system. Indeed, to find the elastic and breakup 
amplitudes one needs only on-shell matrix elements of the transition operator. Each of 
these elements can be found by means of a simple iteration procedure (without complete 
solving the matrix equation (11021)1 with subsequent summation of the iterations via the 
Pade-approximant technique. Note that the same Pade-technique is used in the standard 
approach to the solution of the Faddeev equations [s^ ■ 



Figure 16: The dependence of the GPU-acceleration on dimension of the basis N x N, here N = N, 
for the realistic nd scattering problem with Nijmegen I NN potential at J = dashed line shows the 
acceleration for the step 2 (calculation of the permutation matrix P^), solid line - the acceleration for 
the complete solution (see the main text). 

As was mentioned earlier, the main computational effort in our case is spent on the 
calculation of the free WPs overlap matrix P°. Because these elements are computed 
independently from each other, the algorithm is very appropriate for parallelization and 
implementation on multiprocessor systems, in particular on a graphics processing unit 
(GPU). We have adopted the multi-thread GPU algorithm for a fully realistic calculations 
for nd scattering above the breakup threshold and have attained the great (more than 
tenfold) acceleration in a comparison with the calculation on the same PG without using 
GPU [^. The Fig. [TOl illustrates dependence of the GPU-acceleration in the solution 
of the Faddeev equation with a realistic NN interaction on the dimension of the WP 
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basis. The distinctive feature of the GPU-calculations is that the acceleration grows with 
increasing the basis dimension which is clearly seen in Figure. 


8. 7. Construction of an effective potential for nd scattering 

It has been suggested [32| , that the nd scattering problem can be solved by construct¬ 
ing an effective nd interaction potential. 

Indeed, let us consider again the equation for the transition operator (1961) and divide 
the channel resolvent Gi into two parts (according to the projectors F and Q introduced 
in the Section 7) 

Gi = Gif + Giq. (Ill) 


It is evident that this division corresponds just to bound-continuum and continuum- 
continuum parts of the channel resolvent, i.e. Gif = and Gig = Gf®, where the 
first operator contains channel states including deuteron bound state \zi) . 

Then, instead of a single equation (19611 . one can derive two equations: 


U = Pvi + PviG^^V (112) 

V = Pvi + PviG^’^V. (113) 


It is easy to see that the first equation essentially determines the elastic amplitude and 
is an equation for two-particle scattering, while the second equation defines just the 
effective interaction potential of deuteron as whole with an incident neutron and takes 
into account inelastic processes. 

In the traditional approach, the solution of equation (11131) is practically very hard, 
because it requires knowledge of the channel resolvent in the subspace orthogonal to the 
bound state (i.e. its continuum-continuum part). In the WP approach, the solution of 
this equation presents no difficulty, since the channel resolvent Gf ^ is well known here. 
Then the method for solving the problem of nd scattering by means of effective potentials 
can be represented schematically as follows: 

(i) introduction of the 3WP basis for the channel Hamiltonian and construction of an 
analytical finite-dimensional representation for Gf ® operator; 

(ii) solution of eq. (11131) in the above WP representation; 

(iii) solution of eq. (11121) in a two-body free WP basis using a matrix form of V operator 
obtained at the previous stage. 


9. Summary 

We described in the present paper a general technique for the continuum discretization 
in few-body scattering problem based on projection of scattering operators and respective 
wave functions onto the discrete basis of stationary wave packets. The basic idea behind 
the approach is a similarity (within some restricted space) of exact non-normalized scat¬ 
tering states of the Hamiltonian and the respective discrete and normalized wave-packet 
basis states. So that, such a WP projection allows us to transform the complicated sin¬ 
gular multi-dimensional integral equations describing scattering (like general Lippmann- 
Schwinger or Faddeev-Yakubovsky equations) to regular matrix equations which can be 
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solved directly within computational procedures similar to those used in bound-state 
calculations. 

This novel approach has a few characteristic features which allow to simplify drasti¬ 
cally the whole solution of few- and many-body scattering equations. 

First, owing to some averaging the integral kernels over the momentum cells, all their 
complicated moving singularities are smoothed out and as a result one gets the simple 
matrix equations with finite matrix elements. This allows one to solve the resulting ma¬ 
trix equations directly on real energy axis without any contour deformation, continuation 
to complex energy plane or special interpolation procedures. 

Second, instead of fully off-shell t-matrix at many energies entering the integral kernel 
in a conventional approach, one uses the initial potential (the matrix of which is calculated 
easily) and the matrix of the channel resolvent which is found by means of explicit 
formulas. 

Third, the exact scattering wave packets in pair subsystems are treated as usual 
normalized excited states. This makes possible to consider a three-body breakup process 
as an inelastic scattering to normalized pseudostates. Such a replacement simplifies 
greatly the description of three- and few-body breakup processes. 

Fourth, due to a specific matrix structure of resulting matrix equations one can 
organize the massively-parallel computing with use of the graphics processing unit. This 
makes it possible to carry out all the computations via many thousands of threads on 
the GPU inside a desktop PC and leads to real ultra-fast calculations for the scattering 
problems in few-body systems. 

The developed wave-packet approach to solving few-body scattering problems is uni¬ 
versal and may be used in different branches of nuclear, atomic and chemical physics, in 
quantum statistics and nuclear matter theory. Also it can be directly generalized to a 
relativistic case. 
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Appendix A. Explicit formulas for the channel resolvent eigenvalues 


Appendix A.l. Eigenvalues of two-body resolvent 

Two-body free resolvent eigenvalues for the case of momentum WPs (/(<?) = 1) are 
dehned by the general formula (14711 and thus they have the following form: 


9^{E) 


4(ln 

qd^l q- qi 


+ ln 


g + gi 
g + gj-i 


in [9{q - qk-i) - 9{q - qk)] 




(A.l) 


where q = \/2pE and the combination of the Heaviside 0-functions means that the 
imaginary parts of the eigenvalues don’t vanish only in a single interval to which the 
respective on-shell momentum value belongs: q € 'Sk- 

The energy averaged eigenvalues are defined by integral 

9^{E)Uq, EeDfe. (A.2) 
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Explicit formulas for the g-packet case now are the following: 


5f = 


1 


Dhdi . 


QiV-Qi~^ 


ITT 


(A.3) 


where 


E E {-lf-^'^^-^'[qk-±qA\^\qk'±qA. 


k'—k — l i'—i — l 

In the case of the energy FPs with weight function /(g) = . M the formula for the 


free resolvent eigenvalue is the following: 


while for the averaged one 


E — £i-i 


E-£, 


Di 


[e{E - Ek-i) - e{E - Ek 


k 1 TTA S 


(A.4) 

(A.5) 


where 


Wk^= ^ E - Ee] In \Ek' - Ee 


k'—k—l i'—i — l 


Appendix A.2. Eigenvalues of three-body channel resolvent 

The eigenvalues for the case of energy WPs defined in (ITSh l for the BC part of indices 
are the following 


Re[Gr“'5(A)]fc, = ^ln 


D, 


E-i-i el — E 


■4-E 


lm[G^'^P{E)]k, = -j- {e{E, + e*k - E) - 0{£,., + el - E)} 


(A.6) 


These BC eigenvalues have the same functional form as eigenvalues of the two-body free 
resolvent given in (IA.4I) . The only difference is that bound points £j of partition of the 
free sub-Hamiltonian hj continuum are shifted to eigenvalues of the sub-Hamiltonian 
hi. 

The real parts of the CC-parts of the channel resolvent have the following form : 


I In IA — A_ 


1} 


Re[G^^^iE)]k, = + A-) In |A + A_| + (A - A_)] 

^{(A-l-A+)ln|A-HA+| + (A-A+)ln|A-A+|}, (A.7) 


EkEj 


where: 


A = 


— n 

E*-E, A_=^'= 


A+= 


Dj 


-J - 2 ’ r - 2 

and Dj are energy widths of bins for sub-Hamiltonians hi and /ij respectively. 
The imaginary parts of the CC-eigenvalues also have an explicit analytical form 

TT 


Im[G^“^(£;)]fc, =- 


EkEj 


{(A -F A+)0(A + A+) + (A - A+)0(A - A+) - 

(A + A_)6>(A -h A_) - (A - A_)6>(A - A_)}. (A.8) 
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